Nodal Surfaces in Photoemission from Twisted Bilayer Graphene 
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Selection rules and interference effects in angle resolved photoemission spectra from twisted 
graphene bilayers are studied within a long wavelength theory for the electronic structure. Using a 
generic model for the interlayer coupling, we identify features in the calculated ARPES momentum 
distributions that are controlled by the singularities and topological character of its long wavelength 
spectrum. We distinguish spectral features that are controlled by single-layer singularities in the 
spectrum, their modification by gauge potentials in each layer generated by the interlayer coupling, 
and new energy-dependent interference effects that directly probe the interlayer coherence. The 
results demonstrate how the energy- and polarization- dependence of ARPES spectra can be used 
to characterize the interlayer coupling in twisted bilayer graphenes. 

PACS numbers: 73.22.Pr,78.67.Wj,79.60.-i 



I. INTRODUCTION 

Angle resolved photoemission spectroscopy (ARPES) 
is a well developed tool for mapping out electronic bands 
in solids. Recent applications to single layer and few- 
layer graphenes have demonstrated that, in addition to 
obtaining the energies and momenta of the quasiparti- 
cle bands with high resolution, the measured variations 
of the photoemission intensities also encode informa- 
tion about the phase structure of the wavefunctionglMSl^ 
The matrix elements that couple the initial quasiparticle 
states to the outgoing free electron states are momen- 
tum dependent and produce modulations in the inten- 
sities of ARPES momentum distributions at a constant 
initial state energy. These matrix elements typically con- 
tain a coherent superposition of amplitudes for emission 
from the different sites and layers that are represented in 
the initial state wavef unctions. These in turn depend on 
the (conserved) value of the crystal momentum parallel 
to the surface plane as well as the energy and polarization 
state of the exciting radiation. 

Interpreting these intensity modulations presents a 
complex problem that has been addressed recently for the 
case of single layer and some forms of bilayer graphene^. 
Graphene presents a nearly ideal platform for explor- 
ing this phenomena for two reasons: (i) Graphene is a 
two dimensional material in which there is no disper- 
sion of the quasiparticle bands due to a third (perpen- 
dicular) component of the crystal momentum, (ii) The 
quasiparticle band structure contains point singularities 
around which the internal phases in its wavefunctions un- 
dergo a complete twist. Indeed, interfering amplitudes in 
ARPES from single layer graphene have been observed 
and have striking consequences. The photoemission in- 
tensity on a constant energy momentum space contour 
vanishes along a particular direction (labelled the "dark 
corridor' ; where the emission amplitudes from the two 
sublattices turns out to exactly cancel. The orientations 
of these dark corridors are rotated for photoemission from 
states near symmetry-related zone corner points when re- 



solved in the extended Brillouin zone. For ^B-stacked 
bilayer graphene the emission patterns are more complex 
but they have been analyzed similarly to identify the sign 
of the dominant interlayer Hamiltonian matrix elements 
between neighboring sites in the adjacent layers^. Theory 
predicts that single layer graphene with a strong Rashba 
spin-orbit interactions should have a distinct photoemis- 
sion signature in its spm-resolved ARPES reflecting the 
entanglement of its spin and orbital (pseudospin) degrees 
of freedonP. 

All of these previous theoretical analyses take as their 
starting points a (presumed known) model for the low en- 
ergy electronic Hamiltonian and deduce the consequences 
for the ARPES momentum distributions as a function of 
the initial state energy. In this paper we take the com- 
plementary point of view and consider the ARPES sig- 
natures of a generic Hamiltonian for a graphene bilayer. 
Our approach is motivated by our interest in better un- 
derstanding the electronic structure of "twisted" multi- 
layer graphenes where the symmetry axes of the two lay- 
ers are rotationally misalignect^. It is now widely appre- 
ciated that this misalignment has a profound effect on the 
electronic coherence between neighboring layer^SHMEll^ 
though developing a microscopic theory for it has proven 
to be an elusive goal. Current theories indicate that 
this problem is very rich since kinematical constraints 
define different regions of energy and momentum where 
the two layers can act as "strongly coupled" or "nearly 
decoupled" even for a single structure. These energy- 
momentum sectors depend on the misalignment angle in 
a way that is not yet fully understood. For small angles 
of misalignment theory predicts that new sp ectra l fea- 
tures and narrow bands emerge at low energylSIill, Ex- 
isting experiments are providing conflicting information 
about the nature of the interlayer coupling^ ^^^^^ a dif- 
ficulty that may arise from physical differences between 
twisted graphenes that are made by different experimen- 
tal method^. 

Since both intra-layer and interlayer electronic coher- 



ence in multilayer graphenes can be identified by their 
interference signatures in ARPES, we suggest that deduc- 
ing the coupling from experiment on twisted graphenes 
rather than attempting to calculate it from microscopic 
theory is a promising direction. In this paper we follow 
this approach and present calculations of ARPES intensi- 
ties using a generic model for a twisted graphene bilayer. 
The inputs to the model are the Dirac models for two 
decoupled layers, a momentum offset due to the rota- 
tional misalignment, an electrostatic asymmetry (bias) 
between the layers and a general long wavelength matrix 
that couples the pseudospin amplitudes on neighboring 
layers. From this starting point we calculate the electron 
momentum distributions at fixed energy and the ARPES 
distributions which are weighted by the momentum de- 
pendent photoemission probabilities for various polariza- 
tion states of the exciting radiation. A careful analysis 
of the ARPES momentum distributions shows that the 
"dark corridor" known for single layer graphene changes 
its shape for twisted graphene and allows one to discrim- 
inate between various models for the coupling between 
layers. We present calculations of the ARPES intensi- 
ties that illustrate this effect and a symmetry analysis of 
the relevant matrix that allows us to interpret the spec- 
tra. The results demonstrate how experimental study of 
the ARPES spectra as a function of the energy and po- 
larization of the exciting radiation can be used to fully 
characterize the interlayer coherence in these multilayer 
structures. 

In this paper we first briefly review the long wavelength 
description of twisted bilayer graphene with a theory that 
is appropriate to small rotation angles and in Section II 
we derive an expression for the photoemission transition 
matrix elements based on this model. In section III we 
survey the topological transitions in the band structure 
that occur in this model as the Hamiltonian parameters 
are varied. Section IV shows the calculated ARPES mo- 
mentum distributions showing a rich assortement of in- 
terference phenomena which are unique to twisted multi- 
layer graphenes. Section V gives an analysis of the inter- 
ference effects with in terms of the underlying symmetries 
of the model. We present an analysis of the interference 
patterns in terms of its "energy-flattened" nodal surfaces 
in momentum space which provides a useful diagnostic 
for various forms of interlayer coherence. A further dis- 
cussion of the experimental signatures of these predic- 
tions is given in Section VI. 



low energies there are two electronic bands that touch at 
-E = at discrete points located at the Brillouin zone 
corners vK where v = 1(— 1) denote the K{K') points. 
These bands disperse linearly around the contact points 
and are described by a Hamiltonian at crystal momentum 
k that can be linearized in the difference q = k vK 



Hs = fivpo- ■ {vq^ex + Qyi 



(1) 



where cr are Pauli matrices acting in the space of sublat- 
tice (pseudospin) degree of freedom. 

For a graphene bilayer we adopt a four-component ba- 
sis with a conventional ordering of its four site-layer de- 
grees of freedom (yli, Bi, A2, i?2)- When the symmetry 
axes of the two layers are aligned and the interlayer cou- 
pling is set to zero the Hamiltonian for the two layers is 
simply a doubled version of the single layer Hamiltonian 
and can be written in the form 



Hb = Hs To 



(2) 



where r are Pauli matrices acting on the layer degree of 
freedom. In AB stacked graphene the symmetry axes of 
the two layers are aligned and the two layers are coupled 
through a hybridization of amplitudes on different sub- 
lattice sites in the two layers e.g. (^1, B2) as described 
by the bilayer Hamiltonian 



m 



■crnal 



= 'HsTo 
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i<yxTx — (^yTy) 



(3) 



where 71 is the strength of the interlayer coupling. 



In this work we are concerned with the situation where 
the symmetry axes of the two layers are rotated by a rel- 
ative angle 9. It is convenient to regard each of the layers 
as rotated by angles ±6/2 with respect to a common ref- 
erence structure in which the K point is oriented along 
the X direction as shown in Figure 1. With this conven- 
tion the K points in the two layers shift by momenta 
AK± which are, in complex vector notation 



AK+ = 



(e±»«/2 _ 1) 



K 



(4) 



The components of q are resolved in the original unro- 
tated frame so that the effect of the rotation is to induce 
both a shift of origin and the phases of the complex mo- 
menta appearing in Hb according to the replacement rule 



II. LONG WAVELENGTH MODEL 

A. Small Angle Twisted Bilayers 

Our work uses a long wavelength low energy theory 
of the electronic states in graphene bilayers. In single 
layer graphene the electrons propagate on a honeycomb 
lattice with two sublattice sites labelled A and B. At 



Qx + ky ^ qx,± + ky,± 

<±+<.± - {qx+^qy-{e^''/'-l)K)e^^<^/' 

(5) 

Expressed in a common frame, the offset single-layer 
Dirac operators of Eqn. (IT]) are 



■H± = hvpo- ■ {vq'^j^ex + q'±e 



y,±^v' 



(6) 




qx 



Coefficient Parameterization 


I 


II 


tg (71 - 73)/9 


43.3 


8.3 


CAA 74 + (71 -73)73 


130.0 


69.0 


CAB (71 + 273)73 


130.0 


340.0 



TABLE I: Parameters (meV units) in interlayer hopping 
Hamiltonian for a twisted bilayer fitted to the interlayer in- 
tersite amplitudes 71, 73 and 74 in its Bernal stacked zones. 
Model I: 71 = 390 meV, 73 = 74 = 0. Model II: 71 = 390 
meV, 73 = 315 meV and 74 = 44 meV. 



FIG. 1: The Brillouin zones for each layer (red and blue 
denote different layers) of a twisted bilayer are rotated by 
angles ±^72 with respect to a reference zone (dashed) with 
a corner K aligned with the x axis. In the long wavelength 
model presented in the text all momenta q are measured with 
respect to the reference Brillouin zone corner and resolved 
along the {qx,qy) axes as shown. 



When the rotation angle 9 is small the interlayer coupling 
across a twisted bilayer admits a description analogous to 
the 7 term in Eqn. 3. To see this, before passing to the 
long wavelength limit (Eqn. 1) we partition the bilayer 
lattice Hamiltonian Hiatticc into 2x2 layer diagonal and 
layer off diagonal blocks 



n 



lattice 



T^{r) H2 



(7) 



When 6 is small the registry between layers is modulated 
in a supercell evolving smoothly through local zones with 
{A1B2) stacking, {B1A2) stacking and {A1A2) stacking. 
In this situation T^r) in Eqn. (ItI) is a smooth, periodic 
and local 2x2 matrix function of f acting on the pseu- 
dospins in each layer. The smoothest periodic matrix- 
valued function satisfying these constraints is 



T(rO =io + J2 *" 



^^g,^ 



(8) 



where i„ (n = 0, 6) are 2x2 matrix-valued constants and 
Qn (n = 1,6) are vectors in the first star of superlattice 
reciprocal lattice vectors. The matrix coefficients t„ can 
be further simplified by exploiting lattice and rotational 
symmetries. In a periodic supercell with its AB, BA and 
AA zones centered at positions r^ , r^j and fj respectively, 
the coefficients i„ for the three even elements of Gn are 



p-iG„-r-i g-iG„-r^ 



-iGn-rg 



-iGn-^-y 



tg 



z 1 

z* z 



(9) 



where 



exp(27ri73) and i„odd = t* 



Sublattice 



symmetry requires that the constant matrix ig has the 
form 



to 



caa Cab 
cba cbb 



(10) 



In the small 6 limit, the interlayer coupling is described 
by three real parameters caa = cbb, cab = cba and 
tg. One can choose these parameters by specifying the 
direct (71) and skew (73, 74) interlayer hopping parame- 
ters in the AB registered zones of a twisted bilayer. Ta- 
ble I makes this translation by giving the constants in the 
modulated interlayer hopping model of a twisted bilayer 
in terms of the three interlayer intersite amplitudes. Note 
that if the interlayer hopping amplitudes is long range, 
i.e. if it varies slowly on the scale of the graphene primi- 
tive cell, then tg is small compared to the constant term. 
When the skew hopping terms are isotropic (73 = 74) the 
diagonal and off diagonal terms in the constant matrix 
have the same amplitude {caa — cab)- 

As examples. Table I evaluates these constants for two 
models. In Model I a direct interlayer term 71 is re- 
tained and the skew hopping terms are arbitrarily set to 
zero. Even in this extremely short range model one finds 
that the constant terms dominate the spatially modu- 
lated part of the interlayer Hamiltonian. In Model II 
the constants are fit to the Slonczewski-Weiss-McClure 
(SWMc) parameterization of the interlayer amplitudes in 
Bernal graphite retaining its skew hopping terms^^ and 
as expected the amplitude of the modulated term is even 
weaker. Note also that the SWMc parameterization also 
shows a strong asymmetry between 73 and 74. These 
two amplitudes describe interlayer hopping processes at 
the same range but in different crystallographic direc- 
tions with respect to the graphene symmetry axes. This 
asymmetry, deduced from fitting experimental data, is 
not contained in two-center tight binding theories that 
posit an isotropic interlayer hopping model. The asym- 
metry arises from a crystal field anisotropy in the Wan- 
nier states appropriate to single layer graphene. 

The constant matrix to is the supercell-average of the 
interlayer hopping operatop251. One can distinguish three 
different kinds of behavior depending on whether its off 
diagonal elements are dominant, its diagonal elements are 
dominant or they are equal. The SWMc parameteriza- 
tion realizes the first type of behavior since, as shown in 
Table \n~K\ it is dominated by its site off-diagonal parts: 
Cab ^ CAA- Alternatively, any two center tight bind- 
ing model that assumes an isotropic interlayer hopping 
model is a member of the last family which describes 
a different type of interlayer coherence as we show ex- 



plicitly below. Finally, as noted in our earlier work, the 
absence of a significant Fermi velocity renormalization in 
a family of twisted multilayer graphenes as well as some 
puzzling features in their measured ARPES spectra can 
be understood using the the second class of interlayer 
models in which interlayer hopping on the same sublat- 
tices controls the interlayer Hamiltonian. In this paper 
we will regard the constant matrix to as unknown a priori 
and study the signatures of each of these three forms in 
its ARPES spectra. 

In the four-band description of the twisted bilayer 
problem the interlayer coupling Hamiltonian H+,- = 
(FiCT^: -|-r2cro)ra; where Fi > F2 in the first class, F2 > Fi 
in the second class and Fi = F2 is a marginal state that 
separates them. Allowing for the possibility of different 
electrostatic potentials ±y/2 on the two layers we there- 
fore consider the family of Hamiltonians in the v = 1 
valley 



'Ht = -—- [cr ■ q'_^{To + T^) + cr ■ q'_{To - T^)) 



V 



O-QTz + (FiCTa; + F2Cro)T:E 



(11) 



where the related Hamiltonian for ^' = — 1 is obtained 
by the substitution cr -^ —a*. The Hamiltonian Eqn. 
(Ill is the starting point for the calculations we present 
in this paper. It is a function of the rotation angle 9, the 
electrostatic potential difference V and coupling coeffi- 
cients F„. Choices of F„ allow us to explore the space of 
symmetry allowed interlayer coupling models. 

Note that the projection of the original lattice Hamilto- 
nian H+^- onto the Dirac basis changes the phases of the 
interlayer terms because of the oscillation of the K{K') 
poin t pseudospin basis functions. Thus the interlayer in 
( 11 1 are multiphed by phase factors exp[j(G ■ di — G' ■ d'j)] 

where G{G') are either zero or primitive reciprocal lattice 
vectors in the two separate layers that connect equivalent 
zone corners points and di(d') are the sublattice posi- 
tions in the two layers. Simultaneous boosts by a triad 
of (G, G') pairs transforms the Hamiltonian ( |ll| among 
three pairs of zone corner points in which the twist- 
induced offsets are /S.K and its ±27r/3-rotated counter- 
parts. For definiteness in the remainder of this paper we 
present calculations for the choice G = G' = as given in 
Eqn. (11) describing emission from states near a single 



K point. 



Photoemission Matrix Elements 



In angle resolved photoemission, an incident photon 
with energy hio excites a Bloch electron into an outgoing 
free particle state whose energy Ep = p^ /2me and prop- 
agation direction p = (p||,pz) are measured outside the 
sample. Conservation of energy identifies the initial state 
energy Ei with respect to the work function W 

E,^E„ + W-huj 



and conservation of the parallel component of momentum 
identifies the relative momentum q for emission near val- 
ley vK + G 

q = p\\ vK G 

The incident light is coupled to the electrons through 
an interaction Hamiltonian 



■Hint = ~ev ■ A^ 



(12) 



where the velocity operators are obtained by differenti- 
ating Eqn. (flT|; v^ = dUT/dq^ giving 




2 To T sin 2 T^ 



cos 



± sin I T, 



cos I To 




(13) 



Equation 1 1 3| defines two 4x4 matrix operators acting in 
the pseudospin-layer orbital space for each of the two or- 
thogonal polarizations. These operators are independent 
of momentum q but depend on the rotation angles ±0/2. 
The matrix elements for photoexcitation with photon 
polarization e are 



M{p,q) = {^>\v-e\^g 



(14) 



To evaluate Eqn. ( 14 ) we project the outgoing plane 
wave state f/'p onto the same basis used to represent the 
initial state ipq. Crucially, this projection introduces a 
phase difference 4> = Pzd/h between plane wave ampli- 
tudes on the two layers separated by a vertical distance 
d. For photoemission near the graphene Brillouin zone 
corner, by using excitation energies in the range 30-100 
eV, (f) can be varied over a range of approximately An. 
Note that the parallel component of the momentum pu 
is conserved, laterally phase matching the outgoing state 
to the initial Bloch state. Therefore, under typical ex- 
perimental conditions the analogous lateral interference 
effects are small and nearly energy independent. (They 
would be exactly zero for perpendicular emission.) These 
small effects are not included in our calculations. Thus, 
written in the orbital-layer basis, the final state appear- 
ing in Eqn. ( 14 1 in our model is 



/ 1 \ 



V'^ 



oi</' 



(15) 



\e"*' J 



(151 



The matrix element in Eqn. ( 14 ) is the inner product of 
the operator Eqn. (13) between an occupied eigenstate 
of Eqn. ( 11 ) and the plane wave final state given in Eqn. 



Our analysis presents three different momentum dis- 
tributions for interpretation of the ARPES intensities. 
The first is a map of the spectral function at A{q, Ei) as 
a function of initial state energy Ei unweighted by the 
transition matrix elements 



Aiq,E,) = J2SiEn{q)-E,) 



(16) 



summed over occupied bands n of Eqn. (111. The second 



is a map of a momentum distribution of the simulated 
ARPES intensities I{q,Ei) which includes the state de- 
pendent transition matrix elements 

I{q,E,)=Y, \Mn{p,q)fS{E„iq)-E,) (17) 



As noted above the matrix elements |A^(p, q)P also de- 
pend on the polarization state of the incident light. The 
third is a map that "flattens" the energy resolution of 
Eqn. (17), giving instead for each of the occupied bands 



the distribution of just its momentum- and polarization- 
dependent squared matrix elements 



Pniq) = \Mn{p,q)\' 



(18) 



The summand in Eqn. (17) is the product of the spec- 



tral functions ( 16 1 which expresses the kinematical con- 



straints on the experiment with the transition amplitudes 
( 18 1 that encode the phase information. Maps of the un- 



constrained distributions Pn{q) are useful for exposing 
the quantum geometric origins of the interference pat- 
terns that are accessible in the ARPES intensities. Ex- 
perimentally one can determine the nodal structure of 
( 18 ) by the evolution of a series of spectra that sweep Ei 



through the occupied bands. 



III. TOPOLOGICAL TRANSITIONS 



The spectrum of the Hamiltonian in Eqn. (11) ex- 



hibits topological transitions as a function of its inter- 
layer coupling parameters Fi and r2 as shown in Figure 
2. In the literature on twisted gra phenc s a frequently 
used model treats the case Fi = F j^" * ^^ * ^^ where the long 
wavelength coupling when the interlayer amplitudes are 
isotropic functions of the relative position between two 
sites projected into the xy plane. The spectrum of this 
model shows a pair of Dirac singularities at _E = and a 
saddle point singularity at higher energy where the two 
Dirac cones merge and hybridize. A representative spec- 
trum with these features is shown in the far left panel of 
Figure |2] Perturbation theory in the interlayer coupling 
predicts a ^-dependent renormalization of the Dirac ve- 
locities of its low energy features which are largest in the 
limit of small rotation angles when the two Dirac features 



are nearly congruent. Calculations beyond perturbation 
theory suggest that this low angle regime supports rich 
low energy physics with the emergence of new low energy 
nearly flat bands whose narro w ban dwidth oscillates as 
a function of the rotation angl j^^ l ^ ^. 



Symmetry does not require Fi = F2 and in ear- 
lier worl}2SI we noted that empirical parameterization 
schemes, most notably the Slonczewski-Weiss-McClure 
(SWMc) model suggest that there can be a very strong 
asymmetry between these two parameters. This pro- 
duces striking changes to both the symmetry and the 
topology of the electronic bands both at low energy where 
the Dirac points are effectively decoupled and at higher 
energy where they merge. Using the results of Table 
I, we see that the conventional SWMc parameterization 
describes a situations where Fi ^ F2, i.e. the supercell- 
averaged interlayer coupling is dominated by its sublat- 
tice off diagonal terms. An extreme version of this sit- 
uation one can examine the case Fi 7^ 0,r2 = 0. If 
the coupling strength Fi is weak this model shares many 
features of the isotropic model, with a slightly lower sym- 
metry as shown in the middle panel of Figure [2] In our 
earlier worlPSI ^e found that after a gauge transforma- 
tion this theory can be cast into a form where it describes 
two momentum-offset Dirac cones with opposite helicity 
coupled by a scalar (i.e. sublattice diagonal) interlayer 
coupling matrix. In that analysis one finds that the inter- 
layer coupling generates a gauge potential in each layer 
that displaces the two Dirac nodes towards each other 
as the rotation angle is decreased. Decreasing the rota- 
tion angle increases the effective coupling strength Fi. 
Ultimately, at a critical rotation angle the nodes merge 
and annihilate as shown in Figure [3J For still smaller 
angles two new point singularities emerge describing two 
new Dirac modes of a strongly interlayer-coupled limit. 
It is important to appreciate that the momentum differ- 
ence between the Dirac nodes is not determined purely 
geometrically by the rotation angle, but instead it is gen- 
erally changed by the strength and symmetry of the cou- 
pling between the layers. Also note that in the crossover 
regime the system has a spectrum reminiscent of AB 
stacked graphene, albeit on a reduced energy scale. The 
presence of the singularities at £' = is symmetry pro- 
tected and, as is the case in j4i?-stacked graphene, cannot 
be removed unless the layer symmetry is lifted. 



r 



In the opposite limit (Fi = 0, F2 7^ 0) one finds two 
low energy features that have the same helicity and repel 
each other in momentum space so that no such E — Q 



merger occur regardless of the coupling strength F2 . In- 
stead, at finite energy the two Dirac cones cross to pro- 
duce new symmetry protected points which are second 




FIG. 2: Energy surfaces in the four band model for a twisted graphene bilayer for three different models for the long wavelength 
interlayer coupling matrix: (left panel) Model I with l/int oc (I + ax), (middle panel) (Fi = F2 = 0.05) Model II with Vint oc a^ in 
the weak coupling regime (Fi = 0.13, F2 = 0), (right panel) Model III with Vint oc I (Fi = 0, F2 = 0.15). The blue dots denote 
the positions of the Dirac points of a pair of twisted but decoupled sheets. The actual Dirac points at E — are displaced 
from these points in Models II and III due to gauge potentials generated by the interlayer coupling. 







FIG. 3: Energy surfaces showing a transition from weak to strong coupling in interaction Model II with hvp^K = 0.15. As 
the coupling parameter Fi is increased the renormalized Dirac points evolve from weak coupling (left panel: Fi = 0.13, F2 = 0) 
and merge at a critical point (central panel: (left panel: Fi = 0.15, F2 = 0)) and produce two new Dirac points in the strong 
coupling phase of the model (right panel: (left panel: Fi = 0.20, F2 = 0)). The blue dots denote the Dirac points for a pair of 
twisted by decoupled layers. 



generation Dirac singularities near ±A as shown in the 
right hand panel of Figure [2j Despite the appearance 
of a symmetry-allowed band crossing at the midpoints 
between the E = Q Dirac nodes, the system retains a 
saddle point structure where the density of states is en- 
hanced in an extended region of momentum space that is 
laterally displaced away from the position of the crossing 
point. The low energy features of this type of coupling 
are quite distinct from those of Model I and II. Notably, 
a perturbative velocity renormalization of the low energy 
Dirac nodes is symmetry forbidden for this class of in- 
teraction models. This can provide an explananation for 
the absence of a velocity renormalization inferred from 
Landau level spectroscopjPfil and from ARPES^, both 
of which find a Dirac velocity for twisted BLG that is 
the same as the velocity for single layer graphene within 



experimental error. It also provides a interpretation for 
ARPES experiments that clearly show a crossing of the 
Dirac cones cross instead of a hybridization induced an- 
ticrossing even in the limit of low rotation angles^H. 

The spectra displayed in Figure |2] and [3] show a distinct 
low energy topology controlled by the symmetry of the 
long wavelength interlayer coupling matrix. The singular 
points in these spectra have a striking signature in the 
simulated ARPES momentum distributions as discussed 
below. Additionally there are special lines in momentum 
space where ARPES selection rules are operative and 
these can be used to discriminate between these mod- 
els and even to directly identify the momentum shifts of 
their Dirac feature s due to the gauge coupling from the 
interlayer potential^^'— . The details of this analysis are 
presented in the following sections. 
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FIG. 4: (Top panel) Evolution of the constant energy surfaces for the interaction Model I {hvpAK = 0.15, Fi — .05, F2 — .05) 
taken at energies —0.35, —0.12, —0.07, 0.12,0.35. The far right panel is the result with a interlayer bias. (Middle panel) 
Momentum distributions weighted by the transition probablities for x polarized light. (Bottom panel) Momentum distributions 
weighed by the transition probabilities for y polarized light. 



IV. ARPES MOMENTUM DISTRIBUTIONS 

In this section we discuss the constant energy momen- 
tum distributions calculated for the three limiting models 
discussed in Section III, presented both as unweighted 
distributions and as distributions that are weighted by 
their photoemission transition probabilities. 

The top panel of Figure |4] shows the unweighted mo- 
mentum distribution for the fully symmetric model with 
(Fi = F2 = F) where the interlayer matrix has the form 
Mnt = F(I + (Ti) (Model I). At energies well below the in- 
terlayer scale —hvp^K the constant energy contours are 
self intesecting loops. The inner loop maps out a con- 
stant energy contour on the lowest valence band and as 
the energy is increased (i.e. as it is made less negative) 
this inner loop collapses to a point and vanishes at the ex- 



tremuni of the lowest band. The remaining single surface 
develops a "bean" shape that encircles both Dirac points 
highlighted by the bold dots. At still higher energy this 
remaining contour fissions at energy of the saddle point 
producing two new closed curves that separately encircle 
the two Dirac points. At positive energy these constant 
energy contours expand and reconnect at a positive sad- 
dle point energy to form a single surface. At still higher 
energy these surface "folds" contacting at an energy cor- 
responding to the minimum of its highest energy band 
leading again to a constant energy surface in the form of 
a self intersecting loop. In the presence of an interlayer 
bias (shown in the far right panel) the surfaces main- 
tain the same topology but are inflated/deflated around 
the two Dirac points in response to the layer symmetry 
breaking potential. 



r 



The lower two panels of Figure |4] shows these distribu- 
tions at the same sampling energies but weighted by the 
squared modulus of their transition matrix elements for 
photoemission using light that is x polarized (top row) 
and y polarized (bottom row). For either polarization 
we observe a modulation of the simulated emission in- 
tensities the emergence of a pattern of null points where 
the probability for photoemission vanishes. A symmetry 
analysis of these nodal surfaces is presented in Section 
V. A comparison of these density plots shows that the 
intensity patterns undergo a complex evolution as a func- 
tion of energy with the nodal points undergoing pairwise 
creation and annihilation on constant energy surfaces in 
each band at critical values of the initial state energy. 
As shown in the rightmost two panels, these null points 



are robust features of the simulated ARPES momentum 
distributions and occur even when the layer symmetry is 
broken by an interlayer potential. 

The top panel of Figure [S] shows the unweighted mo- 
mentum distributions using the interlayer coupling model 
V^nt = Ticrj, (Model II). This describes an interlayer 
Hamiltonian that is dominated by its hopping terms that 
connect different sublattices, a situation that is very 
nearly realized by the SWMc parameterization. This 
breaks the symmetry of Model I, so that there are two 
separate (nonintersecting) constant energy surfaces at 
energies below —Hvp^K. As the energy increases the 
innermost loop collapses to a point and vanishes (not 
shown) identifying the band extremum of the lowest en- 
ergy occupied band. At higher energy the remaining 
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FIG. 5: (Top panel) Evolution of the constant energy surfaces for the interaction Model Vint oc a^ {hvpAK — 0.15, Fi = 
0.13, F2 = 0.0) at energies (—0.35, —0.07, 0.12, 0.35) with energy values (left to right). These parameters describe this model in 
the weak coupling limit, i.e. below its topological transition. The far right panel is the result with an interlayer bias. (Middle 
panel) Momentum distributions weighted by the transition probablities for x polarized light. (Bottom panel) Momentum 
distributions weighed by the transition probabilities for y polarized light. 



ring shrinks and fissions at a saddle point energy forming 
two loops which in turn collapse around the two layer- 
coupled Dirac points singularities. As noted in Section 
III these points are generally displaced from the unper- 
turbed Dirac points of the two decoupled layers due to 
a gauge potential generated by the interlayer coupling. 
In the weak coupling limit (i.e. for large rotation angle) 



these new point singularities are located along the ver- 
tical axis of these plots. As the interaction strength is 
increased the system these singular points merge and the 
system undergoes a topological transition to a new strong 
coupling regime where the singular points are symmetri- 
cally positioned along the horizontal axis. 



The lower two panels of Figure [5] show the proba- 
blity weighted distributions calculated for Model II for 
the same representative energies. Again we observe a 
complex modulation of the emission intensities with null 
points identified in the constant energy surfaces for both 
polarizations. The locations of the nodal points in these 



surfaces exhibit an interesting oscillation as a function 
of band index and light polarization. Note for example 
the interchange of the positions of the emission maxima 
and nodal points in the second and third bands under 
exchange of the x and y polarization. 



Figure l6] presents the unweighted momentum distri- 
butions obtained from a similar analysis using a cou- 
pling model of the formViut = r2cro (Model III). In this 
limit the interlayer coupling is dominated by amplitudes 
connecting sites on the same sublattice. A coupling of 
this type would occur in a putative AA stacked geom- 
etry which, though apparently excluded by the conven- 
tional SWMc parameters of Table I, would provide an 
explanation for the absence of a Fermi velocity renor- 
malization and the intersection of Dirac bands observed 



in twisted graphenes grown on SiC (0001). At low en- 
ergy the constant energy surfaces form two disconnected 
loops. The signature of coupling in the form of Model III 
is the collapse of the inner ring at critical energy identi- 
fying second generation Dirac point of Figure [2] and its 
re-emergence as a closed loop immediately above this en- 
ergy (not shown) . As the energy is increased still further 
this ring bifurcates at a pair of saddle points positioned 
displaced along the horizontal (qx) axes. These two sur- 
faces collapse around the two E — Dirac nodes in the 
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FIG. 6: (Top panel) Evolution of the constant energy surfaces for the interaction Model Vint oc I {hvFl\K = 0.15, Fi = 0, F2 = 
0.15) at energy values (—0.35, —0.12, —0.07, 0.12, 0.35). The far right panel is the result with an interlayer bias. (Middle panel) 
Momentum distributions weighted by the transition probablities for x polarized light. (Bottom panel) Momentum distributions 
weighed by the transition probabilities for y polarized light. 



model. Here, because of the gauge potential generated by 
the interlayer coupling the E — Q nodes are repelled from 
the geometrically defined Dirac points. As one crosses to 
positive energies this evolution is reversed leading to the 
reconstruction of two ringlike surfaces for large positive 
energies. 

The lower two panels of Figure l6] show the distribu- 
tions of Model III weighted by their transition probabil- 
ities. Again we observe the complementary positions of 
the emission maxima and nodal points from a single band 
for X and y polarized excitation. 



V. SYMMETRY ANALYSIS 

The nodal patterns identified in Section IV can be un- 
derstood by a consideration of the analytic structure of 
the complex transition amplitudes M.{p, q) in Eqn. ( 14 ). 



The key observation is that in particular circumstances 
identified below this complex transition amplitude degen- 
erates to a real function of q. In this situation M will 
generically vanish along lines in reciprocal space (where 
it changes sign) rather than at points (where it wraps its 
phase). The evolution of the nodal patterns in ARPES 
distributions then correspond to the evolution of constant 
energy surfaces through a network of lines along which 
AA is zero. These lines can terminate at singular points in 
the band structure. In the simplest picture, these points 
can be identified with the rotationally offset Dirac points. 
The actual situation is richer however since, as we have 
seen, these points are generally displaced from the ro- 
tationally defined Dirac points by gauge potentials pro- 
duced by the interlayer coupling and in addition second 
generation symmetry protected singularities can appear 
elsewhere the band structure. The identification of these 



lines, their continuity and termination points in q space 
provides a powerful probe of the nature of the interlayer 
coherence. 



A. Intralayer Interference of the Transition 
Amplitudes 



The layer diagonal blocks in Eqn. 11 are momentum- 
displaced and -rotated versions of the HamiltonianfTland 
have the symmetry 



"Hs = (Jx'H*s(J: 



S"x 



(19) 



so that the eigenfunctions are equal weight states 



x(g) = 



-ia(q) 



\/2 I ±e*"('J) 



(20) 



The velocity operators of Eqn. [13] are layer block- 
diagonal and 0-dependent and can be represented com- 
pactly 



(21) 



in the top layer with positive rotation angle 9 where cr^ is 
the /i-th 2x2 Pauli matrix. Reversing the sign of 9 in this 
expression gives the related velocity operator projected 
on the bottom layer. Combining Eqns. (20 1 and (211 we 
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identify the transition states 



xUq) 



X^iQ) 



xiq) 



x{q) = 



1 



±e'(a(q)-e/2) 

g-j(Q(q)-e/2) 



Te 



i{a(q)-0/2) 



V2 



^-i{a{q)-0/2) 



X — polarized 

y — polarized 

(22) 



The photoemission matrix elements obtained by project- 
ing these transition states onto the plane wave final state 
[15] for X polarized excitation are 

Mt{p,q)^^co4a{q)-e/2) 



mAp.q) 



^/2 



sin(a(q) - 61/2) 



(23) 



and for j/-polarized excitation 

1 



M-{p,q) = ^cosiaiq) - 9/2) 



(24) 



The structure of these matrix elements has been stud- 
ied previously in the context of single layer graphene. 
Note the reversal of the q-space modulations the of x 
and y polarized transition amplitudes for the positive 
and negative energy eigenstates. Importantly the transi- 
tion amplitudes vanish along lines where a{q) = 9/2 or 
a{q) = 9/2 ± tt/2. These nodal lines terminate at points 
where the layer projected Hamiltonian ex I and %"*" and 
X~ are degenerate. Crossing this contact point reverses 
the assignment of ± indices to the upper and lower en- 
ergy bands, thereby terminating the nodal line in a single 
band and transferring a nodal line to its partner. This 
phenomenon is responsible for the one sided "dark corri- 
dor" seen in photoemission from single layer graphene^*^. 
Equations [23] and [24] are rotated versions of that result 
where the new dark corridors are rotated at angles ±9/2 
in the two layers due to the rotational misalignment. 

Any real combination of the transition amplitudes in 
Eqns. (23) and (24) describes linearly polarized excita- 
tion at some general angle of polarization. We see that 
its effect is simply to rotate the nodal line for photoemis- 
sion while preserving an endpoint anchored at a singular 
point of degeneracy. By contrast elliptically or circularly 
polarized light requires a phase lag between two orthog- 
onal polarizations of the radiation. In this situation the 
transition matrix element is necessarily represented by 
a two dimensional parameter space (it must be a com- 
plex scalar function) and it can vanish only at special 
points in momentum space where two independent tran- 
sition amplitudes vanish. Such points would be difficult 
to measure in ARPES since this would require fine tuning 
the initial state to a single critical energy. 



These conclusions hold quite generally for any system 
described by the symmetry of Eqn. ( 19 ) in its layer 
projected Hamiltonians. For coupled BLG, if the in- 
terlayer coupling matrix V respects sublatticc symmetry 
(as occurs in models in the form of Eqn. 



11), one can 



obtain an effective layer projected theory that respects 
sublattice symmetry by integrating out one of the layers. 
Such a system generically exhibits one-sided line nodes in 
the photoemission matrix elements. Since the eigenf unc- 
tions of this layer projected Hamiltonian are equal weight 
states of the form given in Eqn. (20) its matrix elements 



are controlled by the evolution of a single scalar parame- 
ter a{q). This requires that its nodal lines are continuous 
except at point singularities where Va is undefined and a 
nodal line can either terminate or switch between bands 
that are degenerate at the singularity. As noted in Sec- 
tion IV the endpoints of these nodal lines will generally 
deviate from the rotationally defined points of degener- 
acy if a residual gauge potential is induced by integrat- 
ing out the interlayer coupling. Nonetheless the physics 
responsible for this behavior is essentially the singular 
structure of the single layer Hamiltonian perturbatively 
renormalized by interactions between layers. 



B. Interlayer Interference of the Transition 
Amplitudes 

Transition matrix elements derived from the Hamilto- 
nian in Eqn. (11) support additional interference fea- 



tures that arise directly from its interlayer coupling am- 
plitudes. These provide a more sensitive probe of the 
symmetry of the interlayer Hamiltonian in twisted BLG 
and are analyzed in this section. 

Our approach exploits a symmetry in Eqn. ( 11 ) along 
the line Qy = where because of the reversal of the mo- 
mentum offsets zLAK/2 the layer-projected Hamiltoni- 
ans in the top and bottom layers are related 



n 



bottom 



n 



top 



(25) 



This symmetry is broken for a general two dimensional 
momentum q but it remains a local symmetry for all 
momenta qx along the midline. Along this line the eigen- 
functions are equal- weight combinations of the four layer- 
orbital degrees of freedom. Additionally, since the Hamil- 
tonian commutes with the operator R = TxCr^, these 
eigenstates can be indexed by their parities (±) under 
R and take the form 



/ 



-ZQ( 



vl/±(g) 



.ta{q) 



e 

^giaiq) 



(26) 



The four-component velocity operators are 
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/l 



Va 



ie/2 



-ie/2 



CTp 

<7., 



r 



/I 



.-ie/2 



je/2 



(27) 



and the transition states are obtained by application of 
these operators to the initial states VP*, giving for x po- 
larized light: 



/ 



*^ 



V^-^^ 



\ 



and for y polarized light 



V 



Vy^^ 



gi(a-e/2) 

p-i(a-9/2) 

_^g-i(a-e/2) 

V ±e'("-«/2) y 



_g-i(a-e/2) 
_^g-i(a-e/2) 



(28) 



(29) 



The matrix elements for photoemission in this case de- 
pends on the relative phase cj) for the layer amplitudes in 
the outgoing state in Eqn. (15). For a; polarized light 
the matrix elements are 



M + {p, q) ^ cos(a - 61/2) cos(0/2) 
M^{p, q) = -i cos{a - 9/2) sin(0/2) 

and for y polarized light 



(30) 



M + {p, q) = sin(a - 6/2) sin(0/2) 
My{p,q)^is\n{a-e/2)cos{(l)/2) (31) 



In Eqns. OO^ and (31) the phases a and 9 are indepen- 



dent variables (the former is determined by coefhcients in 
the Hamiltonian and the latter is the rotation angle) so 
that generically a— 9/2 is not a multiple of 7r/2. Thus the 
relevant interference physics is fully controlled by the in- 
terlayer phase (p. When the outgoing waves from the two 
layers are in phase (</> = 2?Ti7r) emission from states with 
even i?-parity is allowed for x polarized and forbidden 
for y polarized light. Conversely, emission from the odd 
i?-parity states is allowed for y polarized light and for- 
bidden for X polarization. When the phase lag between 
layers (p = {2m -f l)7r these selection rules are exactly 
reversed. As is the case for intralayer interference, nodal 
lines are absent for circular or elliptical polarization since 
this requires a simultaneous vanishing of the matrix ele- 
ments in two orthogonal excitation channels. 

When ^ (0, tt) the emission between layers is neither 
exactly in or out of phase and nodal lines do not occur 
along the midline. Instead they are replaced by troughs 
(local minima) in the ARPES momentum distributions. 



The (^-derivatives of the momentum distributions mea- 
sured for two orthogonal linear polarizations allow one to 
extract the internal phases a that define the wavefunc- 
tions ^± along this line. Physically these derivatives 
can be measured by measuring the differential change of 
the emission intensity at each position in q space dis- 
tribution as a function of the excitation energy. Using 
Eqns. (30) and (31) one finds that the intensities for 



even i?-parity state in two orthogonal polarizations are 



h = cos^{a-9/2) 
ly = sin2(a-6l/2) 



1 + cos(0) 



1 — cos(0) 



(32) 



and therefore the ratio ryx 

(energy) independent and allows one to identify a{q) 



-{dly/d(j))/(dh/d(j)) is 



a{q) = — h arctan 



\Av(9) 



(33) 



Eqn. (33) holds everywhere along the midline when 
the interlayer bias V — 0. When the interlayer bias is 
nonzero a similar expression can be used to determine the 
wavefunctions along a hyperbolic locus in q-space where 
the wavefunction ^± is an equal weight state shared be- 
tween the two layers. Illustrative examples are given in 
Section V.D. 

Eqns. (30) and (31) predict that nodal lines at gj, = 
extend over the entire midline when the R- parity of the 
initial state is constant along this line. This is always the 
case when the bands are nonintersecting and illustrations 
of this are given in the following sections. Note however 
that when a symmetry protected band crossing occurs 
along the midline, the i?-parity changes its sign at a sin- 
gular point of degeneracy identifiable by the termination 
a nodal line in one band (ordered by its energy) and 
its appearance in another (also ordered by energy). The 
Qy = midline is perpendicular to the offset AK between 
the bare Dirac nodes of the two layers so that crossings 
of this type do not occur in the simplest model of the 
rotated bilayer (Model I). However, in refined models for 
twisted bilayer graphene symmetry protected crossings 
do occur, notably along the line Qy — near E — in 
the strong coupling limit of Model II and at the energy 
of the second generation Dirac point (near the interlayer 
scale hvpAK) for any coupling strength in Model III. 
This provides a unique spectroscopic diagnostic of singu- 
larities in their band structures. 
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The interlayer selection rules can be understood from 
the transformation of the current operators V^ under the 
symmetry operation R. Ke(j,) have even(odd) i?-parities, 
and for cj) — 0(7r) the outgoing plane wave state in Eqn. 
(15) is even(odd). Thus the even parity initial state ^+ 
can be excited into the outgoing final state with ip — 0{tt) 
only using x{y) polarization. Similar considerations ap- 
ply for general interlayer phase difference (/) with a rota- 
tion of the principal axes as noted above. 



C. Branching Nodal Surfaces 



The nodal patterns discussed in Section V A arise from 
cancellation of photoemission amplitudes from the de- 
grees of freedom in the individual layers, and in Sec- 
tion VB| they arise from cancellation of amplitudes from 
different layers. There are special points in q space 
where both destructive interference conditions are sat- 
isfied. These can be identified as branching points where 
a nodal line running along the midline qy — devel- 
ops branches that continue to the (renormalized) Dirac 
points of the two layers. The locations of these branch- 
ing points can be determined from consideration of the 
bilayer coherent wavefunctions along the high symmetry 
q-y = midline and are analyzed in this section. 

Using the Hamiltonian in Eqn. (11) we partition the 



Hamiltonian at qy ^ into a Hamiltonian along the mid- 
line at qy = and a qy dependent piece 



niy) 



nqy 
nqy 



0)+H'{qy) 
0) + qyToa2 



(34) 



Therefore the eigenstates near the midline can be ex- 
panded in terms of the eigenstates along the midline 



^ 7n \Qy 



*m(g.y =0) 

ly}^ *n(0) 7^ ^ (35) 



E„ 



En 



The midline eigenfunctions are also eigenfunctions of i? = 
<^xTxj and since H' is odd under R the sum connects 
only states of opposite i?-parity. In a similar way the 
photoemission matrix element for the m-th band in the 
/x-th polarization Mfj, can be expanded 



Mjn,t,iqy) = A^m,M(0) 



qy 



E 



(^^bMl*«(0)>(*n(0)|(T2To|*,„(0)) 



(36) 



The sum in Eqn. ( 36 ) runs over the two states that 



reverse the R parity of the initial state ^m(O) and using 
Eqns. (30) and ( [3l| ) when Aim,ii{0) = one sees that 
Ain,ti (OjTs nonzeroror each of these admixed states. This 



means that the photoemission matrix element turns on 
linearly as a function of the transverse momentum qy 
and the nodal surface is an unbranched line for a general 
momentum q^. 

However a singular point can occur along the midline 
{qb, 0) where the two contributions in the sum cancel. To 
locate such a point we label the two energy denominators 
between the state m and its two intermediate states ni 
and 712 by Ai = E^n — En-^ and A2 = -£„— i?„.2- Then the 
Qy-linear changes to the matrix element can be written 
as a real function of the phase angle a giving 



5M 



qy 



cos 2a cos(a — 6/2) sin 2a sin(a — 9/2) 



Ai 



= BqyCOs{l3-a + e/2) 



(37) 



where /3 — arctan((Ai/A2) tan2a) and i? is a constant. 
A critical branching point occurs when the right hand 
side of Eqn. ([37]) is zero. A nodal line forms two branches 
above and below the qx axis at such a critical point, once 
formed these arms continue to the Dirac points along the 
qy axis. This can be understood as a consequence of the 
continuity of the intralayer nodal lines that must connect 
the two Dirac points on opposite sides of the q^ axis. 



D. Examples 

1. Model I 

The nodal lines for Model I with 14nt oc (I -I- cr^;) are 
shown in Figure [7] We observe that the nodal lines that 
occur along the midline qy = propagate without ter- 
mination. This occurs because the presence or absence 
of a nodal line is determined by the polarization state of 
the light and by the i?-parity of the state, which does 
not change as a function of qx along the midline. The 
interlayer nodal lines also occur as complementary pairs, 
i.e. they are present for a particular band in x polariza- 
tion only when absent in y polarization and vice versa. 
This occurs whenever the phase lag (f> for emission be- 
tween layers is a multiple of tt. The phase difference (f> 
can be tuned continuously by varying the excitation en- 
ergy and for 7^ mn the nodal lines are absent in both 
X and y polarizations. However for intermediate (/> these 
lines are recovered when the principal polarization axes 
are also rotated according to the formulas given in Eqns. 
(30) and (31 ). No such energy dependence occurs for the 



nodal surfaces that terminate on the Dirac points which 
arise essentially from intralayer interference in the tran- 
sition amplitudes and are therefore i/i- independent. The 
energy and polarization dependence of the nodal lines can 
therefore be used as a smoking gun to discriminate be- 
tween intralayer and interlayer interference phenomena 
in the momentum distributions. 
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FIG. 7: Nodal lines in the photoemission amplitudes calculated for interlayer coupling Model I with Vint oc (I + ax)- The top 
row is for x-polarized light and the bottom row is for y polarized light. The four columns are for the four bands of the model 
ordered by energy. The x polarized amplitudes in bands 2 and 3 have one-sided nodal lines at constant values of qy while the 
y polarized amplitudes have nodal curves due to a gauge potential generated by the interlayer coupling. In bands 2 and 3 the 
nodal lines terminate at points (red dots) that can be identified with the Dirac points of the decoupled layers (green dots). 



There are additional nodal lines in the second and third 
bands that terminate at the Dirac points. The termina- 
tion and the appearance of these features correspond to 
a transferring a line node between the bands that touch 
at the Dirac points. The signature of the nodal patterns 
produced by the coupling in Model I is the appearance, 
in X polarized excitation, of one-sided horizontal nodal 
lines that terminate at the Dirac points. This occurs 
because the interlayer operator in Model I has one zero 
eigenvalue, implying that there is one state with a partic- 
ular pseudospin polarization in each layer that cannot be 
transported to its neighboring layer. This state is the an- 
tisymmetric eigenfunction of a^ and it cannot be coupled 
to the outgoing cTa; -symmetric final state by x polarized 
radiation. Since these states are localized to the individ- 
ual layers, they are eigenfunctions of the single layer 2x2 
Dirac operators for which the pseudospin polarization re- 
mains constant along radially directed lines in q space. 
Thus there are two such states with nodal lines that ter- 
minates at the unperturbed Dirac nodes of the two layers. 
In y polarization the related lines nodes share the same 
termination points but the pseudospin polarization is a 
function of both g^, and Qy and is constant along a curve 
in momentum space, as shown in the lower panel. The 



curvature of this line is a manifestation of the momen- 
tum dependence of the gauge potential produced by the 
interaction with the neighboring layer. The termination 
points for these lines are the same for x and y polarized 
excitation, and thus the point singularities in this model 
can be identified with the geometrically determined Dirac 
points of two decoupled layers. This is a special feature 
of this class of interaction models. 
2. Model II 



In Model II V^nt oc ax so that the interlayer coupling is 
controlled by terms that connect opposite sublattice sites 
on neighboring layers. The signature of this model is the 
existence of a weak and strong coupling limit separated 
by a critical state where the renormalized Dirac points of 
the two layers are merged. The nodal structure in these 
two regimes is illutrated in Figures [8] and |9] The dimen- 
sionless coupling strength Fi = cab/^vf^K so that the 
strong coupling limit can be realized for sufficiently small 
rotation angles. Using the data of Table I, one expects 
the strong coupling regime to describe the physics for 
rotation angles 9 < A° . 



The nodal structure for Model II in the weak coupling 
regime is shown in Figure [8J It shares some features in 
common with Model I and the nodal pattern in y po- 
larization is quite similar. The main difference is in x 
polarization where "flat" one-sided nodal lines of Model 



I evolve into dispersive nodal curves in Model II. As noted 
above the curvature of these lines is a momentum-space 
manifestation of the interlayer gauge potential for the 
twisted bilayer which is absent by symmetry for a single 
polarization in Model I. 
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FIG. 8: Nodal lines in the photoemission amplitudes calculated for interlayer coupling Model II with Vint oc ax- The top 
row is for x-polarized light and the bottom row is for y polarized light. The four columns are for the four bands of the model 
ordered by energy. The data are shown for Model II in its weak coupling regime before the merger and reconstruction of the 
E = Dirac points. The renormalized Dirac points (red dots) are displaced from the geometrically defined Dirac points of 
two decoupled layers (open green circles). The nodal lines in bands 2 and 3 terminate on the renormalized Dirac points. An 
additional nodal line along the midline is produced by interlayer interference of the transition amplitudes. 
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FIG. 9: Nodal lines in the photoemission amplitudes calculated for interlayer coupling Model II with Vint oc a^ in the strong 
coupling regime. The top row is for x-polarized light and the bottom row is for y polarized light. The four columns are for the 
four band of the model ordered by energy. The data show the effect of reconstruction of the E = Dirac points (red dots) in 
the strong coupling regime which are positioned along the ±qa: axis with the geometrically defined Dirac points (open green 
circles) along the zbg^ axis. The nodal lines in bands 2 and 3 terminate on these renormalized Dirac points. An additional 
nodal line along the midline occurs due to interlayer interference of the transition amplitudes. 



In the strong coupling regime the nodal patterns in 
Model II shown in Figure|9]are more interesting. Here the 
reconstructed strong coupling Dirac points occur along 
the high symmetry g^. axis. Thus the interlayer nodal 
lines can and do terminate at singular points along the qy 
axis. Interestingly these spectra also support branching 
points as shown. The nodal lines that branch away from 
these branch points propagate without termination to 
large momenta as shown. 

3. Model III 



connect the same sublattice in the two layers. The nodal 
patterns calculated for this model are shown in Figure 
[TOl The signature of this family of interaction models is 
the appearance of second generation symmetry-protected 
Dirac points at finite energy. 



In Model III the interlayer coupling matrix I^nt oc I 
so that the coupling is dominated by amplitudes that 
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FIG. 10: Nodal lines in the photoemission amplitudes calculated for interlayer coupling Model III with Vint oc I. The top row 
is for x-polarized light and the bottom row is for y polarized light. The four columns are for the four bands of the model ordered 
by energy. Red dots denote the renormalized E — Dirac points in the model that are displaced away from each other along 
the ±qy axis by gauge potentials generated by the interlayer coupling. The open green circles denote the undisplaced Dirac 
nodes of two decoupled layers. The nodal lines in bands 2 and 3 terminate on these renormalized Dirac points. An additional 
one-sided nodal line terminating at a second generation Dirac point occurs along the midline occurs in all four bands due to 
interlayer interference of the transition amplitudes. 



Interlayer nodal lines propagating along the Qx axis 
collide with these second generation points where they 
switch between bands. The fundamental signature of this 
class of models is the appearance of new one-sided nodal 
lines along the midline qy = 0. This behavior is clearly 
evident in both polarizations in Figure |10[ The nodal 
lines that terminate on the renormalized E — Dirac 
points are again curves rather than straight rays, mani- 
festing the q-dependence of the gauge potential generated 
by the interlayer coupling. 



4- Phases along the rmdlme 

Figure ITT] shows the internal phases a in the wave- 
functions ^± calculated along the midline qy = using 
Eqn. (33 I for the three interaction models discussed in 



the previous sections. We have verified that the the inter- 
nal phase a determined by differentiation of the emission 
intensities is independent of the value of the phase lag 
(j) for emission from the two layers. A careful inspection 
of these plots allows one to further discriminate between 
the various interaction models. 

All the models are characterized by an evolution from 
a = to a = 7r/2 as a function of increasing q^. The 



phase angle a is half the phase difference between sub- 
lattice amplitudes in the same layer. Thus this evolu- 
tion describes the transition from states composed of the 
asymptotic (large \qx\) layer eigenstates (1,1) to (1,-1) 
which is a common feature in the negative energy solu- 
tions in all three interaction models. 

The intermediate behavior is quite different in these 
models and can be used to analyze the symmetry of the 
interlayer coupling in the intermediate and small qx re- 
gion. For example, in Model II (middle panel) the inter- 
layer coupling matrix V^nt oc ax and using Eqn. (25 1 the 



phase a is then identified as the geometrical Dirac angle 
arctan(Aiir/g2,). Note that this is the same for both neg- 
ative energy bands at all values of qx in this model. By 
contrast Model III (right panel) shows a different phase 
angle a for the two negative energy states as a conse- 
quence of their interlayer coupling. The signature of this 
family of interaction models is both this phase splitting 
and a jump discontinuity that occurs at qx = where the 
two negative energy bands touch and the i?-parity of a 
given band changes. Model I combines features of both: 
the phase angle a deviates from the Dirac angle defined 
by the rotational misorientation and the angles a evolve 
smoothly from to 7r/2 in both manifolds. 



r 



VI. DISCUSSION 

In this paper we have seen how the energy and polar- 
ization dependence of ARPES momentum distributions 



from bilayer graphene can be used to measure the phase 
structure of its layer-coupled wavefunctions. The nodal 
surfaces discussed in this paper can be measured in two 
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FIG. 11: Variation of the internal phase parameter a as a function of q:c along the midline qy = 0. The three panels are for 
three limiting coupling models: left panel: Vint en I + (Tj,, middle panel: Vint oc a^, right panel: Vint oc I. The solid(dashed) 
curves give the variation of the phase angle a in the first(second) occupied band. 



ways. One may sweep the excitation energy and keep 
the outgoing electron energy fixed, in which case one 
measures the momentum distribution as a function of 
a varying initial state energy. This protocol effectively 
sweeps the constant initial state energy contours through 
a pattern of nodal lines in their transition probabilities 
and allows a determination of momentum space nodal 
structures shown in Figures 7-10. The evolution of the 
momentum distributions during this type of sweep are 



illustrated in Figure 12 for interaction Models II and III. 
Alternatively, one can sweep the excitation energy and 
detection energy simultaneously, thereby measuring the 
momentum distribution for a single initial state energy. 
This protocol probes the symmetry of single initial state 
wavef unctions. Both methods are useful. The finger- 
print of intralayer interference are nodal patterns which 
do not change as a function of the excitation energy, and 
this seems best suited to the first method. By contrast, 
interlayer interference effects are energy dependent and 
our analysis exploits the differential energy dependence of 
the emission intensities in orthogonal polarizations. This 
requires following the emission from the same intial state 
as the excitation energy is varied and can be analyzed 
most easily by the second method. 

It will be useful to first confirm these predicted en- 
ergy dependences experimentally. We expect there will 
be a clean separation between intra- and inter- layer in- 
terference effects seen in the experimental spectra with 
linearly polarizated light and a suppression of these ef- 



fects for circularly polarized light. Then focusing on the 
linearly polarized spectra, observation of the curvature of 
the intra-layer nodal surfaces that terminate on the pri- 
mary Dirac points will provide a definitive experimental 
signature of the momentum dependence of the gauge po- 
tentials produced by the interlayer coupling. Conversely, 
observation of 0-periodic variations of the emission inten- 
sities along particular lines in momentum space identify 
those features that are controlled by the interlayer coher- 
ence of the electronic states and differential dependence 
of their intensities on the excitation energy can be used to 
determined the sublattice symmetry of the matrix valued 
interlayer coupling potential. 

Although we have focused on the situation with no ver- 
tical electrostatic potential difference between the layers 
the results can be generalized to graphene bilayers with 
scalar layer asymmetry. In this case the relevant momen- 
tum space contours deform to a family of mutually or- 
thogonal confocal hyperbolas and ellipses. The extension 
of our methods to that situation is relatively straightfor- 
ward and will likely be needed for a quantitative analysis 
of experimental data. The methods developed in this pa- 
per are also quite general and can be implemented with 
more sophisticated models for the electronic structure of 
BLG. It is hoped that measurement and analysis along 
these lines will provide a clean experimental resolution 
of the nature of the electronic states in twisted bilayer 
graphenes. 
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